home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / TMTA.CPP < prev    next >
C/C++ Source or Header  |  1995-01-11  |  3KB  |  98 lines
  1.  
  2. //#define WANT_STREAM
  3.  
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmatap.h"
  8.  
  9.  
  10. /**************************** test program ******************************/
  11.  
  12. void Print(const Matrix& X);
  13. void Print(const UpperTriangularMatrix& X);
  14. void Print(const DiagonalMatrix& X);
  15. void Print(const SymmetricMatrix& X);
  16. void Print(const LowerTriangularMatrix& X);
  17.  
  18. void Clean(Matrix&, Real);
  19.  
  20.  
  21. static void process(const GeneralMatrix& A,
  22.    const ColumnVector& X1, const ColumnVector& X2)
  23. {
  24.       Matrix B = A;
  25.       LinearEquationSolver L(A);
  26.       Matrix Y(4,2);
  27.       Y.Column(1) << L.i() * X1; Y.Column(2) << L.i() * X2;
  28.       Matrix Z(4,2); Z.Column(1) << X1; Z.Column(2) << X2;
  29.       Z = B * Y - Z; Clean(Z,0.00000001); Print(Z);
  30. }
  31.  
  32.  
  33.  
  34. void trymata()
  35. {
  36. //   cout << "\nTenth test of Matrix package\n";
  37.    Tracer et("Tenth test of Matrix package");
  38.    Exception::PrintTrace(TRUE);
  39.    int i; int j;
  40.    UpperTriangularMatrix U(8);
  41.    for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
  42.    Matrix X(8,6);
  43.    for (i=1;i<=8;i++) for (j=1;j<=6;j++) X(i,j)=i*j+1.0;
  44.    Matrix Y = U.i()*X; Matrix MU=U;
  45.    Y=Y-MU.i()*X; Clean(Y,0.00000001); Print(Y);
  46.    Y = U.t().i()*X; Y=Y-MU.t().i()*X; Clean(Y,0.00000001); Print(Y);
  47.    UpperTriangularMatrix UX(8);
  48.    for (i=1;i<=8;i++) for (j=i;j<=8;j++) UX(i,j)=i+j+1;
  49.    UX(4,4)=0; UX(4,5)=0;
  50.    UpperTriangularMatrix UY = U.i() * UX;
  51.    { X=UX; MU=U; Y=UY-MU.i()*X; Clean(Y,0.000000001); Print(Y); }
  52.    LowerTriangularMatrix LY = U.t().i() * UX.t();
  53.    { Y=LY-MU.i().t()*X.t(); Clean(Y,0.000000001); Print(Y); }
  54.    DiagonalMatrix D(8); for (i=1;i<=8;i++) D(i,i)=i+1;
  55.    { X=D.i()*MU; }
  56.    { UY=U; UY=D.i()*UY; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
  57.    { UY=D.i()*U; Y=UY-X; Clean(Y,0.00000001); Print(Y); }
  58. //   X=MU.t();
  59. //   LY=D.i()*U.t(); Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
  60. //   LowerTriangularMatrix L=U.t();
  61. //   LY=D.i()*L; Y=D*LY-X; Clean(Y,0.00000001); Print(Y);
  62.    U.ReDimension(8);
  63.    for (i=1;i<=8;i++) for (j=i;j<=8;j++) U(i,j)=i+j*j+5;
  64.    MU = U;
  65.    MU = U.i() - MU.i(); Clean(MU,0.00000001); Print(MU);
  66.    MU = U.t().i() - U.i().t(); Clean(MU,0.00000001); Print(MU);
  67.  
  68.    // test LINEQ
  69.    {
  70.       ColumnVector X1(4), X2(4);
  71.       X1(1)=1; X1(2)=2; X1(3)=3; X1(4)=4;
  72.       X2(1)=1; X2(2)=10; X2(3)=100; X2(4)=1000;
  73.  
  74.  
  75.       Matrix A(4,4);
  76.       A(1,1)=1; A(1,2)=3; A(1,3)=0; A(1,4)=0;
  77.       A(2,1)=3; A(2,2)=2; A(2,3)=5; A(2,4)=0;
  78.       A(3,1)=0; A(3,2)=5; A(3,3)=4; A(3,4)=1;
  79.       A(4,1)=0; A(4,2)=0; A(4,3)=1; A(4,4)=3;
  80.       process(A,X1,X2);
  81.  
  82.       BandMatrix B(4,1,1);  B.Inject(A);
  83.       process(B,X1,X2);
  84.  
  85.       UpperTriangularMatrix U(4);
  86.       U(1,1)=1; U(1,2)=2; U(1,3)=3; U(1,4)=4;
  87.         U(2,2)=8; U(2,3)=7; U(2,4)=6;
  88.               U(3,3)=2; U(3,4)=7;
  89.                     U(4,4)=1;
  90.       process(U,X1,X2);
  91.  
  92.       LowerTriangularMatrix L = U.t();
  93.       process(L,X1,X2);
  94.    }
  95.  
  96. //   cout << "\nEnd of tenth test\n";
  97. }
  98.